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We apply a large-deviation method to study the diffusive trajectories of the quadrature operators 
of light within a reservoir connected to dissipative quantum systems. We formulate the study of 
quadrature trajectories in terms of characteristic operators and show that in the long time limit 
the statistics of such trajectories obey a large-deviation principle. We take our motivation from 
homodyne detection schemes which allow the statistics of quadrature operator of the light field to 
be measured. We illustrate our approach with four examples of increasing complexity: a driven 
two-level system, a 'blinking' three-level system, a pair of weakly- coupled two-level driven systems, 
and the micromaser. We discuss how quadrature operators can serve as alternative order parameters 
for the classification of dynamical phases, which is particularly useful in cases where the statistics 
of quantum jumps cannot distinguish between such phases. The formalism we introduce also allows 
us to analyse the properties of the light emitted by quantum jump trajectories which fluctuate far 
from the typical dynamics. 
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I. INTRODUCTION 

This paper addresses questions about the nature of dy- 
namical phase transitions and crossovers in open quan- 
tum systems from the perspective of observers able to 
make measurements on the environment. We take our 
motivation from experiments using homodyne detection 
schemes p]-[3[ and study the time series or trajectories 
of quadratures of the light emitted from open quantum 
systems. We go beyond recent work [4] which focused on 
trajectories of quantum jumps analyzed from the point 
of view of the counting statistics of, for example, photons 
entering the environment. In this approach, identifying 
the average rate of photon emission as a dynamical order 
parameter allowed to uncover dynamical phase transi- 
tions and crossovers in a number of systems [El-Hof. A 
central idea in the works upon which we build was the 
introduction of a field conjugate to the number of emitted 
quanta, the so-called "s-field" . In equilibrium statistical 
mechanics, the observation of phase transitions is depen- 
dent on choice of the external field used to tune across a 
transition and the use of an appropriate order parameter; 
in this work we demonstrate that different dynamical or- 
der parameters and conjugate s-fields allow exploration 
of different dynamical phase transitions. Specifically we 
show that crossovers between dynamical phases may be 
observed from quadrature measurements which could not 
be found by counting photons, and vice versa. We also 
extend our study to consider the effects of two s-fields, so 
that we can influence both light quadratures and quan- 
tum jumps. This enables us to understand the relation 
between inferences from different dynamical observables; 
for example, we can find the typical quadrature measure- 
ments which would be made in rare quantum-jump tra- 
jectories, when either fewer or more photons are emitted 
than average. 

Our approach takes as a starting point the basic ob- 
servation that (real-time) dynamics is more than stat- 



ics. Sometimes dynamical behaviour, including transi- 
tions and crossovers between different dynamical regimes, 
can be understood from the properties of the stationary 
state. But a more common occurrence is that dynami- 
cal fluctuations are only revealed via (often high-order) 
time-correlation functions, and not by static or one-time 
observables. This suggests that for a proper statistical 
analysis of the dynamics of open systems a 'statistical 
mechanics of trajectories' is needed. Such an approach is 
the so-called s-ensemble method, which has proven use- 
ful for the study of classical many-body systems display- 
ing complex cooperative dynamics such as glasses [1 ll — 
Il5j . The aim is to describe dynamical phases in terms of 
strictly dynamical order parameters and to classify these 
dynamical phases and the changes between them using 
a mathematical and conceptual framework analogous to 
that of equilibrium statistical mechanics. This ensemble 
method for dynamics can therefore be thought of as a 
'thermodynamics of trajectories' approach [il TlTL [l6l - [l8| . 

Dynamical phase transitions are not limited to classical 
systems; transitions and crossovers have been discovered 
in a number of driven open quantum systems. Famous 
examples include the laser [19] , whose behaviour close to 
threshold resembles a closed thermodynamic system in 
the neighbourhood of a continuous phase transition, and 
the micromaser [20|, Hl| , where a flux of atoms is coupled 
to an optical cavity mode and drives the cavity through 
a series of crossovers. However more recently, dynamical 
phase transitions have also been discovered in a broad 
range of fields, spanning decohering spin channels [22j , 
information transport in complex systems (23[, current 
fluctuations in isolated diffusive systems Q and even the 
decoherence of tunnelling molecules [25| . The s-ensemble 
method was recently applied to the study of ensembles of 
quantum-jump trajectories in open quantum systems 
Alt hough s imilar in spirit to ideas in full counting statis- 
tics |26l-[29| , the s-ensemble approaches the problem from 
a different perspective and has now been developed for a 
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FIG. 1. A schematic diagram of the general Markovian system 
plus environment models we study. The interaction between 
the system and the environment is of form cb^ + bc\ where 
the operators (b) are creation (annihilation) operators of an 
harmonic oscillator bath and the operators c and are system 
operators which are related to the Lindblad projectors of the 
system. 



variety of problems where dynamical phases are classified 
according to the counting statistics of quantum jumps 0- 

03,123. 

In this paper, we formulate the application of the s- 
ensemble to light quadratures by means of reduced char- 
acteristic operators. Measuring quadratures of emitted 
light allows the quantum state of light emitted from 
an open system to be probed in a way not possible by 
counting photons. For example, it provides an under- 
standing of whether the light is in a coherent state, a 
squeezed state or another more exotic state. We ap- 
ply this method to a selection of open quantum sys- 
tems, ranging from few- level quantum-optical systems to 
the micromaser, all of which have a stucture illustrated 
schematically in Fig. HJ In contrast to counting the num- 
ber of emitted photons, as in studies of quantum-jump 
trajectories, we consider the temporal accumulation of 
quadratures of light emitted from a system into its envi- 
ronment. We consider this a quadrature trajectory. Un- 
like the quantum-jump trajectories, the accumlation of 
light quadratures in the environment is a diffusive pro- 
cess. However we can define a quadrature activity, de- 
fined as a time-averaged light quadrature, and use this 
as a dynamical order parameter to characterise dynami- 
cal phases in the space of quadrature trajectories. 

Our motivation for studying the statistics of quadra- 
tures is inspired by the experimental technique of homo- 
dyne detection. The X-quadrature trajectories of emit- 
ted light are directly related to the homodyne current [H- 
@], as described in Fig. [2j (It is also sensible for us to 
study other quadratures, since these are accessible via a 
change in the driving Hamiltonian, which we will discuss 
later in Sec. IIII1 ) The statistics of quadrature trajecto- 
ries provide a more natural probe of the system dynamics 
than quantum-jump trajectories and we will show that 
such measurement schemes allow exploration of dynam- 
ical phases in a variety of systems. 

Beyond examining dynamical phases identified by 
quadrature activity, we construct marginal probability 
distributions for general quadrature operators at all an- 
gles in phase space. Using these marginals, we recon- 
struct Wigner distributions Q via the inverse radon 
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FIG. 2. Shown (left) is a setup for homodyne detection and 
(right) a typical photocurrent, /homo, plotted against time 
in such a measurement. A simple homodyne scheme con- 
sists of an input beam and an output beam which is split 
between a detector and a strong coherent field which acts 
as a local oscillator. If we denote &out as the output lower- 
ing operator and S as the local oscillator strength, the total 
transmitted field is bo ut = 8 + frout- Looking at the pho- 
tocurrent I p — dbQ Ut bo ut /dt one may define the homodyne 
current /homo = — £ in the limit of infinite local oscillator 
strength. By expanding I p one finds that the homodyne cur- 
rent is a direct measure of the quantity b-\-b\ which is directly 
related to the X-quadrature. Ref. 0] contains further details 
about the various homodyne detection schemes which allow 
X-quadrature measurements. In this work, we use the input 
and output fields to form the reservoir we depict in Fig. [1] 



transform, to find the state of the emitted light. We 
extend our studies to examine the typical quadrature tra- 
jectories of systems biased by the number of photon emis- 
sions. This allows us to understand the nature of light 
emitted from a system for all quantum-jump trajecto- 
ries, whether rare or typical with respect to the number 
of emitted photons. More generally we study the tra- 
jectories of a particular observable after having biased 
the system towards rare trajectories of another (generally 
non-commuting) observable. This technique also allows 
us to identify the appropriate dynamical order parame- 
ters to understand the dynamical phases in different open 
systems. 

The paper is structured as follows. In the following 
section we introduce the formalism of the s-ensemble 
and describe the Ito calculus methods which we em- 
ploy for the stochastic description of the environment 
we use in this paper. We further develop generalised 
master equations for the dynamics of quadrature-biased 
systems and systems biased towards rare trajectories of 
both quantum-jumps and quadratures. We also discuss 
use of the Wigner function to analyse our results. In sub- 
sequent sections we present our results. In Sec. IIIII we 
discuss a simple driven two- level system and in Sec. HYl 
we extend our analysis to a driven three-level system. In 
Sec. El we study a pair of coupled two- level systems and 
demonstrate that measuring quadrature trajectories al- 
lows identification of dynamical phases which are hidden 
in photon-counting experiments. Finally, in Sec. [VU we 
give a concise account of quadrature trajectories in the 
micromaser, providing insights via a mean- field theory as 
well as exact numerical diagonalisation. In Sec. IVIII we 
give our conclusions. 
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II. FORMALISM AND THEORETICAL 
BACKGROUND 

Several recent studies have applied the s-ensemble 
method to study the statistics of quantum trajectories in 
open quantum systems using a thermodynamic descrip- 
tion [3l], [32j . In this section we will give an account of the 
s-ensemble formalism for quantum trajectories. We dis- 
cuss the specific examples of counting photons and mea- 
suring the quadratures of light entering the environment 
around an open quantum system. First, in Section Hi Al 
we discuss the general application of the s-ensemble in 
the study of the statistics of an observable Q. We then 
reformulate these principles so their application to open 
systems is clear. Section Hi Bl discusses the open quantum 
system methods which we use, which are applied to the 
s-ensemble in Sections III CI and III Dl 



A. An s-ensemble for open quantum systems 

We begin by considering a particular measurable quan- 
tity Q associated with a quantum trajectory of a quan- 
tum system coupled to a reservoir. (For example, Q could 
be the number of emitted photons in a time t, or other 
quantities which we define below.) Projecting the system 
density matrix p(t) on to the subspace where Q takes a 
particular value, one may define a reduced density ma- 
trix p[^ . The probability of such a realisation occuring 
in time t is then given by P t (Q) = Tr [p iQ) (t)] and, after 
long times when the system reaches a steady state, this 
takes a large-deviation (LD) form: 



Pt(Q) 



-t<KQ/t) 



(1) 



where (j)(Q/t) contains all the information about the 
probability distribution of Q at long times. Alternatively 
the statistics may be described introducing a moment 
generating function associated with these probabilities. 
This method proceeds by introducing another density 
matrix p s (t) defined by the Laplace transform 



p s (t) = J p( Q \t)e- sQ dQ (2) 

from which we define the moment generating function 

Z t (s) = Tvp s (t) = J P t {Q)e- s ®dQ. (3) 

In both ([2]) and (|3]) , the integrals should be replaced with 
sums if Q is a discrete quantity, as is the case when count- 
ing photons. A LD form 



Z t (s) ~ e 
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(4) 



— mm s (0Q(s) + (Q/t)s). The full statistics of Q are there- 
fore contained within Oq(s). The density matrix p s (t) de- 
scribes, when s ^ 0, rare trajectories where Q is far from 
the mean. In counting processes, where Q is bounded 
from below by zero, the rare trajectories can be sepa- 
rated into more active trajectories when s < and less 
active trajectories when s > 0. These correspond respec- 
tively to trajectories with a larger or a smaller number 
of counts than in the s = physical dynamics. 

There are two immediate advantages to this s-ensemble 
approach. Firstly, as the LD function 0q(s) is the largest 
real eigenvalue of a superoperator which generates the s- 
biased dynamics of p s (t) (which we derive in the following 
sections), evaluation of the statistics of trajectories can 
be straightforward. Of particular interest are the first 
and second moments, q s and Aq 2 , which may be found 
directly from the LD function via 
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is found at long times, with LD functions 0q(s) 
and (j)(Q/t) related by Legendre transform, <fr(Q/t) = 



where the s subscripts indicated that the expectation val- 
ues are taken with respect to the ensemble of trajecto- 
ries biased by e~ s ®. Secondly, this method provides a 
thermodynamic formalism for non-equilibrium processes. 
The LD functions 0q(s) and (j)(Q/t) are analogues of free 
energy and entropy densities, with s the conjugate in- 
tensive field to the time- extensive Q. Furthermore, the 
so-called activity q s may be used as a dynamical order pa- 
rameter to distinguish dynamical phases, whose bound- 
aries may be crossed by tuning the parameter s, or other 
system parameters. Indeed, such phase boundaries are 
crossed when an s-derivative of 0q(s) becomes discontin- 
uous. 

So far, the s-ensemble studies of open quantum sys- 
tems have explored the thermodynamics of quantum- 
jump trajectories associated with counting quanta emit- 
ted from a system into a Markovian environment. In 
these cases the quantity Q is, for example, the number 
of photons emitted (i.e. the number of quantum jumps), 
K, from a system in unit time. Defining reservoir lad- 
der operators b and in the Heisenberg representation, 
the quantity K corresponds integrating over time the ob- 
servable b^b . Biasing these trajectories, we obtain a LD 
function 6k( s ) from which we extract the dynamical or- 
der parameter k s = (K) 3 jt. 

In this paper, we develop the s-ensemble further to 
study the statistics of the quadratures of light entering a 
Markovian bath from various quantum systems. Choos- 
ing Q to be the X- and F-quadratures of the light, corre- 
sponding to bath operators (6 + &^)/2 and i(b — tf)/2, we 
will use corresponding quadrature activities x s = (X) s /t 
and y s = (Y) s /t when applying s-fields to bias the dy- 
namics towards rare quadrature trajectories. In much 
the same way as the photon activity k s can be associated 
with the average number of photons emitted in a time t, 
the quadrature activities x s and y s are time averages of 



the X- and F-quadratures of the light leaving the system. 
We will gain yet further insight into the dynamics of the 
systems we study by examining general quadratures of 
the form 

X a = cosa X + sina Y (7) 

where the angle a in phase space is illustrated in Fig. [3j 
For these general quadratures, we can consider marginal 
distributions P t (X a ) associated with general quadra- 
tures, which can be extracted from the relevant LD 
functions 6x<*(s) and their s-derivatives using Eqs. (|5j) 
and (j6j). 
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FIG. 3. The quadrature operators define axes of an optical 
phase space. The generic quadrature operator I a may be 
viewed as a rotation of the X-quadrature axis as shown above. 

In the following subsections we will introduce the 
mathematical formalism which allows us to construct s- 
biased ensembles and show how the appropriate LD func- 
tions are derived. We will also explore doubly-biased 
ensembles where we study the typical trajectories of one 
observable for ensembles of trajectories biased by another 
observable. We develop this formulation in Section Hi D[ 
but we first discuss how the LD functions for (singly-) 
biased ensembles are derived. 



are described by Ito increments which obey a quantum 
Ito calculus [3|- We next sketch our method within this 
formalism; we give full details of the derivation in Ap- 
pendix [XJ 

For the rest of the paper we consider the reservoir to 
be an unsqueezed vacuum [33|. We introduce temporal 
increments dB^ (t) and dB (t) which are defined in terms 
of the reservoir ladder operators b(t) and b^(t) by 

/t+dt pt 
b\t')dt' - / b\t')dt', 

/t+dt nt 
b(t')dt'- / b(t')dt'. (9) 

and are evaluated in the immediate future. The dynamics 
are found using a stochastic density matrix Rd) whose 
Ito increment is related to the system Hamiltonian, H, 
the reservoir increments defined above and the Lindblad 
operators of the system (Li): 

dp (t) = -i[H, p]dt - X - Y^{L\Lu P }dt (10) 

+J2 dBl (*) L ^(*) L \ dB (*) 
(*) u P (t) 

+ J2p(t)dB(t)L\. 

i 

Tracing out the reservoir degrees of freedom, one obtains 
the familiar Lindblad master equation 0, IHj for the 
(trace-preserving) system dynamics: 



B. Generalized Master equations, Ito Calculus and 
Quadratures 

Previous studies on quantum-jump trajectories 0, [1| 
have demonstrated that one may identify the LD function 
0q(s) as the largest eigenvalue of an s-modified master 
equation obeyed by p s , 

P s = W S (P% (8) 

where W s is the generator for the s-biased dynamics. We 
will outline how to construct formally these generalized 
master equations in terms of an Ito calculus and char- 
acteristic operators, before proceeding to apply them to 
the study of quadratures in Sec. Ill CI 

We consider the system to be weakly coupled to a 
reservoir in the Markovian regime [3]; within this pic- 
ture the reservoir is treated as a bath of harmonic os- 
cillators which effectively act as a white noise source to 
our system. This type of reservoir admits a stochastic de- 
scription whereby the effects of reservoir ladder operators 
b(t) and 6+ (t), with commutator [b(t),tf(t')] = S(t - t') 



= -i[ff sP0 ] - y,{4^p } + E L ^° L l ( n ) 

i i 

where {•,•} is the anti-commutator and p° = TrR es (p) 
is the density matrix of the system, with TrR es denoting 
the trace over the reservoir Hilbert space. 

To modify this trace-preserving scheme to the desired 
s-biased generalized master equation for trajectories of 
Q, we introduce a characteristic operator, Vq(t). The 
corresponding stochastic increment dQ is related to the 
reservoir increments in Eq. (j9]). The characteristic oper- 
ator is defined by 

V^(t)=exp(-sJ dQit')^. (12) 

From this operator we identify s-biased reduced density 
matrix as 

p s =Tr Res (V$(t)p(t)) . (13) 
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Using Eq. ([T3]) one may derive Eq. (j8j) by finding the 
increment of p s in terms of the increments of Vq and 
p. Then, application of the Ito calculus of dQ in terms 
of reservoir increments (|9j) allows the trace to be taken 
over the reservoir degrees of freedom (see Appendix |A|) . 
This method of deriving W s relies upon the existence 
of a complementary stochastic reservoir operator to the 
system operator of interest. We discuss these operators 
next for the cases of quantum jumps and quadratures. 

We turn our attention to statistics of quantum-jump 
trajectories, before moving on to quadrature trajectories, 
where Q = K. A quantum-jump trajectory is the time 
record of the number of jump events K over a time t. We 
formulate the stochastic process associated with jump 
trajectories in terms of an increment dK. In terms of 
the reservoir increments this jump increment is defined 
as dK = dB^dB/dt which has eigenvalues of and 1 
within the interval [t,t + dt]. Although dK tells us how 
many photons were emitted into the bath in the inter- 
val [£, t + dt], it tells us nothing else about the form of 
the light emitted. We address questions about the quan- 
tum state of emitted light by looking at the statistics of 
the quadrature trajectories of the light emitted into the 
reservoir. We use the X- and F-quadratures to define 
coordinate axes of the optical phase space (see Fig. [3]). 
We will also study the general quadrature operator where 
the observable is Q = Xf , defined in Eq. (0, where a 
is the polar angle. The stochastic increments associated 
with these operators, dX a (t), are defined by 

dX a (t) = i (e- ia dB(t) + e ia dB\t)) . (14) 

Evaluation at a = and tt/2 yields the increment for 
the familiar X- and F-quadratures. Although dX a (t) 
takes the same form for all a, the form of the stochas- 
tic process dK associated with the jump trajectories is 
very different. Therefore, for the rest of this paper we 
use notation where K is biased by a conjugate field s', 
while the quadratures X a are biased by s. With this no- 
tation, we define characteristic operators of the form (fT2j) 
associated with each ensemble of trajectories, 

V£ a (t)=exp(-s f dX a {t')^j 

V^{t)=e X p(-8jdK{f)y (15) 

With these definitions we will construct generalized 
master equations for each ensemble of trajectories in 
Sec. Ill C[ and we will proceed further in Sec. Ill Dl to ask 
about the typical statistics of one process in the biased 
ensemble of the other. 

C. An s-Ensemble for Quadrature Trajectories 

We now turn to the statistics of the trajectories asso- 
ciated with the quadrature operators. We consider pro- 
jections of the density matrix p(t) on to the subspace 



where X a takes a specific realisation. By this we mean 
that the time-integrated quadrature for light entering the 
bath up to time t takes a specific value X a . We de- 
note a corresponding projected reduced density matrix 
by p\ X ^ . As described in Section III Al this is related 
to an s-biased reduced density matrix p s (t), defined in 
Eq. (j2j) with Q = X a , via the Laplace transform 

p s (t) = J p( xa \t)e- sXa dX a . (16) 

which upon taking the trace gives us the moment generat- 
ing function associated with these diffusive probabilities 
Z t (s) = Tr p 3 (t) ~ e t6xa ( s \ where the LD form is valid at 
long times. The activity associated with the quadratures 
is x a = (X a )/t. This, and the second moment Ax a2 are 
found from the derivatives of the LD function Ox™ (s) de- 
fined in Eqs. © and © with Q = X a . This LD function 
is identified as the largest real eigenvalue [3l|, [36| of the 
generalized master equation 

p s (t) = W s (p s ) 

= C (p s ) -E| (e~ ia L iP s + e ia p s L]) + ^p* . (17) 

i 

We derive this generalised master equation using 
Eqs. ([TO]) and (fT2j) . where tracing out the environment is 
performed within the Ito calculus formalism, the details 
of which are presented in Appendix [A] When 5^0, the 
superoperator collapses to the trace-preserving Liouvil- 
lian £, Eq. ([TT]) . Away from s = the s-field biases the 
dynamics towards rare trajectories of the system. The 
properties of s-modified master equations are discussed 
in Ref. [1, Q . Before presenting our results for various 
model systems, we discuss doubly-biased ensembles. 



D. Quadratures of Quantum Jump Trajectories: 
Doubly-Biased Ensembles 

Thus far we have provided the theoretical formalism 
for an s-ensemble of quadrature trajectories. Before pre- 
senting results for a variety of systems, we introduce a 
further interesting problem in this section concerning the 
trajectories of already-biased ensembles. Consider a sys- 
tem biased towards rare trajectories with, for example, 
more quantum jumps than the s f = average. We now 
ask what are the quadratures for these trajectories which 
are, for example, more active with respect to the num- 
ber of emitted photons? In order to answer this question 
we need to introduce two counting fields, s' conjugate to 
the number of emitted photons, K, and s conjugate to 
quadratures X a . 

The derivation of doubly-biased ensembles is an ex- 
tension of the derivation in Sec. Ill C( using Eq. ([TP]) 
with characteristic operators and tracing out the reser- 
voir. The characteristic operators associated with the 
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quadrature increment and photon increment are respec- 
tively V Xa it) and V K (t). With these we can define a 
new s and s' biased density operator, p ss , which incor- 
porates information on the statistics of both ensembles 
of trajectories, 



P ss (t) = Tr r 



v K v Xa v K 



Pit) 



(18) 



The ordering of the characteristic operators is important 
due to the non-commut ability of the two observables. As 
defined, Eq. (fT8j) is interpreted as biasing the jump tra- 
jectories and within these s'-biased ensembles measuring 
the quadrature realisations. We show in Appendix lAl 
how Ito calculus is used to calculate the Ito increment of 
this new double-biased density operator. (We formally 
trace out the reservoir using the Ito tables in Eqs. (|A2|) 
and (|Aip .) We find this new doubly-biased master den- 
sity operator obeys a new generalised master equation 
which reads 



p ss \t) = C (V s ') + (e~ s ' - l) kc P ss 'J 



(19) 



S\/Ke 2 



(e- iot cp ss ' + e i( *p ss 'J 



Notice that if s, s' — >• then the left-hand side reduces 
to the (trace-preserving) Liouvillian. Once again we may 
encapsulate the total long time trajectory statistics us- 
ing the largest eigenvalue, 0k,x<* (s, of the generator 
of these statistics. Furthermore we note that 6k, x a (s, s f ) 
necessarily reduces to the LD function for s'-biased jump 
trajectories, 9k(s'), in the limit s —> [1,0]. To examine 
the typical quadratures of quantum-jump biased systems 
we examine the derivatives of 6k,x<*(SiS')i with respect 
to 8, evaluated in the limit s — >• 0. Note that similar 
doubly-biased ensembles may be created by choosing dif- 
ferent characteristic operators in Eq. ([T8]) . For example, 
switching the characteristic operators in Eq. ([15]) allows 
us to find the statistics of quantum jumps in quadrature- 
biased ensembles. 

Next, we put the formalism into practice in the fol- 
lowing sections with studies of few-level quantum-optical 
systems and a system with more complex dynamics, the 
micromaser. 



III. DRIVEN TWO-LEVEL OPEN SYSTEMS 

In this section we present our results for a simple 
quantum-optical system consisting of a dissipative two- 
level system driven by a laser, shown schematically in 
Fig. HJ We begin by examining the X- and F-quadrature 
statistics of the system and compare it with s-ensemble 
studies for the jump trajectories studied in |4]. Af- 
ter examining both analytic and numerical forms of the 
LD function we turn our attention to the doubly-biased 
statistics of the system. We finish by constructing plots 



of time-independent probability distributions, at various 
s', associated with the X a which we refer to as phase 
portraits before examining the Wigner functions of this 
system. 
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FIG. 4. Schematic diagram of a laser-driven two-level system 
coupled to a vacuum reservoir. 

For the driven two- level system, the generalized master 
operator, W s , is of the form in Eq. (fTTj) . where now we 
have one pair of Lindblad operators L and L\ where L = 
y^c, with k the decay rate and c the lowering operator 
|0) We drive the system with a laser polarised in the 
x-direction, which introduces a driving term ex <j x in the 
two-level system Hamiltonian, 



H = ft (c - 



(20) 



where Q is the Rabi frequency. Here we consider the 
specific choice n = 4ft as this has been shown to be an 
interesting parameter choice in previous works [4], and 
for further reasons we discuss below. We study the X- 
and F-quadratures specifically. In this regime although 
the LD function does not take on a straightforward form, 
writing W s in matrix form we numerically diagonalise it 
and identify its largest real eigenvalue as the LD func- 
tion. Although no concise analytical form exists for the 
LD function, we may Taylor expand the LD functions 
about s = to see how the statistics of the physical 
system behave. For X- and F-quadratures we find the 
expansions 
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(21) 



Firstly we examine the form of the LD functions close 
to s = 0. From the Taylor expansion (|2T|) . we see that 
that $x (s) to O (s 9 ), is symmetric about s = 0. Imme- 
diately we find that the X-quadrature activity for phys- 
ical 5 = dynamics is £ s =o = 0- This symmetry is not 
present in the F-quadrature statistics as the Hamilto- 
nian breaks the symmetry between X and Y\ the odd 
powers of s in the Taylor expansion (f2T|) imply that near 
s = 0, 0y{s) is asymmetric such that, close to the s = 
physical dynamics, the sign of s is relevant. Evaluating 
the physical F-quadrature activity we find y s =o = —2/3, 
the absolute value of which corresponds to the value of 
the physical photon emission rate, k s '=o, as shown in 



FIG. 5. A plot of the F-quadrature activity \y s | and minus the 
quantum- jump activity — k s > for the two- level system when 
biased by the appropriate fields s and s' . Note that the order 
parameters coincide for the physical dynamics at s — s' — 0. 



Fig. [5l This result is specific to the system parameters 
here. While the choice n = AVt was identified as special in 
Ref. [4] due to a self similarity of quantum-jump statis- 
tics, we find this is also the choice of system parameters 
where the F-quadrature and quantum-jump activities co- 
incide. That this relation is found for F-quadratures in- 
stead of X-quadratures is entirely due to the chosen laser 
polarisation. 

Switching our focus to the full numerical forms of the 
LD functions, we show in Figs. [Gal and I6bl that the min- 
imum of 6y is shifted away from s = 0. However, while 
0x(s) and 0y(s) are distinct near s = 0, at very large 
Ox and By both tend towards s 2 /S. While the laser 
breaks the symmetry between X- and F-quadratures, the 
rare trajectories become equivalent again at large \s\. We 
can understand this as the s-field dominating the non- 
equilibrium dynamics over the self-Hamiltonian of the 
system. 

We now explore how biasing the photon activities with 
a counting field s' affects the s = measured quadra- 
tures. The dynamics of our doubly biased system are 
described by the superoperator W ss ' defined in Eq. ([T9]) . 
We diagonalise this superoperator to find the LD func- 
tion 0k,x<x(s, and evaluate the quadrature activities 
by taking the derivative with respect to s at s = 0. 
We look far into both the photon inactive and active 
regimes discussed in Sec. Ill Al by choosing photon count- 
ing fields s f = 5 and —5 respectively. To aid visualisa- 
tion of how the biasing affects the light, we construct 
a time-independent phase portrait of the light and find 
the corresponding Wigner distribution, which we define 
in Appendix [B] To construct our phase portrait we de- 
termine the LD function for equally spaced values of a 
in the range [0,7r] and use the Legendre transform to cal- 
culate e~^ x \ from Eq. (pQ). We plot these probability 
distributions in an activity phase space using coordinates 
x s =o and y s =o (hereafter we drop the s = subscript). 
These plots, or portraits, illustrate the set of marginal 
distributions for all a: the marginal probability distri- 
butions correspond to cuts through the origin at angle 
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FIG. 6. Plots of the LD functions and associated activities 
and variances for the laser-driven two-level system. Shown in 
(a), the LD function for the X-quadrature statistics is com- 
pletely symmetric about s = 0. It is also worth noting that 
as \s\ becomes larger, the dynamical variance tends to 1/4, a 
consequence of the quadratic form of Ox at large \s\. In (b) we 
show that the y-quadrature LD function is now asymmetric 
close to s = 0, while at large |s| it adopts the same quadratic 
form as 0x(s). Both LD functions exhibit a positive activity 
for s < and a negative activity for s > 0. 



a to the abscissa. Our results are shown in Fig. [7J Us- 
ing the values e~^ x ) at each s' ', we also construct the 
Wigner distribution, shown in Fig. [8j with a numerical 
implementation of the inverse Wigner transform. 

Figure [7] shows that the physical dynamics have a 
"heart-shaped" portrait centred on (x,y) = (0,-2/3). 
Although it may appear as though there are bi-modal 
regions, all e~^ x ^ are in fact non-Gaussian unimodal 
distributions where, as a increases from to 7r/2, the 
mean shifts from to —2/3. Interestingly the photon in- 
active system appears more like a vacuum as it is centred 
closer to the origin at (x, y) ^ (0, —0.289). However, al- 
though this is a photon-inactive region, there is always a 
finite probability of quantum jump provided s is finite. 
In a complementary fashion the s' = — 5 portrait "blows 
out" as the centre of each e~^ x ) moves rapidly away 
from zero, when a increases from to n/2. Correspond- 
ing features are shown in the Wigner functions plotted in 
Fig. [8l The Wigner function becomes more concentrated 
about the origin as we enter the photon-number inactive 
region, confirming the notion that the output field be- 
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FIG. 7. Phase-space portraits displaying the marginal probability distributions e~^ x ^ for the two-level system at various 
photon biases s' . We demonstrate that, as we make the system more photon-active, the plot moves away from the origin in the 
negative ^/-direction. In (a), s' — —5 and the center of the plot is (0, -1.534); in (b), s' — and the plot center is (0, -0.667); in 
(c), s — +5 and the plot center is at (0, -0.289). Note that although the axes are labelled x and y, they are not time extensive. 
This is formally equivalent to considering activities when t — 1. 



(b) s'= 



(c) s'= +5 



4 

3 
2 

> 

"o 

2-1 

-2 
-3 
-4 



1 — i i^^i "i — r 



J I I I I L 



35 
30 
25 



-4-3-2-101234 
X-Activity 



20 f 1 
15 
10 



o 

5-1 



-2 h 
-3 
-4 



~i — i — i — i — i — i — r 



o 



j i i i i i i_ 



-4-3-2-101234 
X-Activity 



0.007 
0.006 
0.005 
0.004 
0.003 
0.002 
0.001 


-0.001 



4 
3 
2 

£ 1 

> 

"o 

5-1 

-2 
-3 
-4 



"i — i — i — i — i — i — r 



o 



J i i i i i i_ 



-4-3-2-101234 
X-Activity 



0.0014 

0.0012 

0.001 

0.0008 

0.0006 

0.0004 

0.0002 



-0.0002 



FIG. 8. Reconstructed Wigner functions for the same system parameters as in Fig. a is incremented by O.Olrad over the 
range [0, tt] in the reconstruction. The values s' = —5, and +5 are plotted in (a), (b) and (c) respectively. Focus on the central 
regions where the Wigner function is large: the regions far from the central maxima which take negative values are necessarily 
sensitive to the choice of a-increment used. We have reproduced the Wigner functions with a broad range of a- increments and 
see no changes to the central regions when choosing increments between O.Olrad and O.OOlrad. 



comes more vacuum-like as we proceed to more positive 
s f . 



phase portraits of the system. 



IV. DRIVEN THREE-LEVEL OPEN SYSTEMS 

We now turn our attention to a more complicated 
quantum optical system, the three-level system, de- 
picted in Fig. [9] We drive this system by two resonant 
lasers with Rabi frequencies Qi and ^2 between the |0) 
level and upper levels |1) and |2) respectively. When 
Qi ^> 0,2 the photon emission trajectories are intermit- 
tent 0, [37], [HI; the |1) — >• |0) transition is the active 
or light line transition while the |2) is an inactive level. 
Previous studies |4[ found that the k s > exhibits a dynam- 
ical crossover at s' = 0. The crossover is from an active 
phase s' < 0, where the dynamics are dominated by the 
|1) and |0) levels, to an inactive phase s f > 0, where the 
system dynamics are dominated by long periods where 
the 1 2) level is occupied. With this in mind we will ex- 
amine firstly the quadrature trajectories, then proceed 
to examine the effects of double biasing and examine the 




FIG. 9. A three- level system coupled to a vacuum reservoir 
and driven by two resonant lasers with Rabi frequencies Qi 
and Q2. When Qi ^> Q2 the photon emission trajectories are 
intermittent; the |1) — >• |0) transition is the active or light line 
transition while |2) is an inactive level. 



Considering the three-level system driven by two 
lasers, the generalized superoperator W s is given in 
Eq. (fTT|h still with one set of Lindblad operators L and 
]J . In this system, L = |0) (1| and there are no other 
Lindblad operators because of the null decay rate be- 
tween 1 2) and |0), as illustrated in Fig. [9] However our 
Hamiltonian for this system is slightly different: 
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2 

ff = $^fi<(ci + ct), (22) 

2 = 1 

where a = |0) (i| and cj = |i) (0|. When fti > ^ 2 typ- 
ical photon emissions are intermittent displaying both 
"bright" and "dark" periods (37|, [38|. The LD function 
(s) for each a is obtained by direct diagonalization of 
W s and we examine specifically the case n = AVt i and 

We first study the ensembles of trajectories with X- 
quadrature biases, shown in Fig. HDTa). Here, Ox (s) is 
symmetric in s and exhibits no sharp crossovers. On 
the other hand, the ensemble of F-quadratures exhibits a 
sharp peak in its dynamical variance as we cross s = 0, as 
shown in Fig. fTOT b). This peak corresponds to a crossover 
between two distinct dynamical phases corresponding to 
the changeover of the behaviour of y(s) in this region. 
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FIG. 10. (a) The LD function, activities and dynamical vari- 
ances for the X-quadrature statistics of the three- level system. 
This is completely symmetric about s = 0, but it is not the 
same as the two- level LD function except in the large \s\ limit. 
This can be seen in the dynamical variance, which tends to 
1/4 in this limit, a consequence of the quadratic form of Ox 
at large \s\. (b) The corresponding plots for y-quadratures 
for the three-level system are shown. The activity exhibits a 
marked crossover around s = 0, which appears as a rounded 
step. At the crossover, the corresponding dynamical variance 
is shown to have a sharp peak. 

In the photon emission study in Ref. [4| , the crossover 



in k' s was attributed to a change in the effective behaviour 
where the active side where s' < was argued to be sim- 
ilar to a two-level system. In the inactive phase found 
for s f > 0, the |2) level is occupied for long time periods 
such that few photon emissions occur. We attribute the 
crossover observed in y s to the change in dynamics as- 
sociated with the crossover in photon emission discussed 
in Ref. [4j. To support this claim we examine phase- 
space portraits of the system and the effects of the s'-bias 
on these portraits, by solving the doubly-biased master 
equation (fT7|) . 

First we look at the photon-inactive dynamics with 
s' = —5. The phase portraits and Wigner functions are 
shown Fig. [TTT a): these appear identical to the two-level 
system biased by the same s'-field, shown in Fig. [7^a). 
Moving through to s f = +5, shown in Fig. ITTTc), the 
phase portrait appears almost like a vacuum. This is 
not a true vacuum; both photon activity and y s for this 
portrait are small but non-zero. However the forms of the 
portrait and Wigner function do suggest that for large s', 
the 1 2) state is occupied for a large fraction of the time. 
We illustrate our finding more clearly in Figs ll2al and ll2b[ 
where we show how a contour plot of the s = two- 
level portrait with three-level shifted s = +5 contours 
closely resembles the three- level system s = plot. This 
demonstrates that the physical dynamics of the three- 
level system is effectively made up of an active two-level 
system plus an inactive two-level system (this is where 
the dynamics are dominated by occupations of the |2) 
state) with some "mixing" of the two. 

To determine whether or not it is this k s > crossover 
which is responsible for the crossover in y s we examine 
the statistics of typical jump trajectories for rare quadra- 
ture trajectories. We repeat the procedure outlined in 
Sec. Ill Dl but now we first bias quadratures and then 
count quantum jumps by introducing the doubly-biased 
density matrix 

p*'\t) = Tr res (V s W 8 'V**p) . (23) 
Then, we generate a new master equation 

p s ' s (t) = C (/ s ) + (e~ s ' - l) Kd/'cl (24) 

„2 SyJH \ e~ s ' + 1 ) , N 
+ j/ 8 4 {e- ia ci/ s + e Vet) . 

Solving for the LD function, 0x<*,k{s' \ s) for both the X- 
and Y- quadratures we find k s '=o for various quadrature 
biases s. The typical activity k shows no sharp features 
as we bias the X-quadrature, irrespective of the sign of 
s;. However, if we bias the F-quadrature, the activity k 
and variance Ak 2 show photon-number inactivity when 
s < 0. That is to say, more positive y s leads to smaller 
k. This is demonstrated in Fig. [13l This supports our 
previous assertion that the activities y s and k s > may be 
used as equivalent order parameters in this system. For 
completeness, we show in Fig. [2] how y s =o varies as a 
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FIG. 11. Phase-space portraits of the marginal probability distributions e~^ x ^ at various photon biases s for the three- level 
system. As we make the system more photon-active the plot moves away from the origin in the negative ^/-direction. In (a), 
s' — —5 and the center of the plot is (0, -1.534); in (b), s' — and the plot center is (0, -0.4) and in (c), s' — +5 and the plot 
center is at (0,0). 




FIG. 12. (a) A contour plot of the, s = 0, two-level system 
density plot plus the three-level system, s' = +5, (effectively 
a one-level system due to inactivity) density plot with its cen- 
ter shifted by 0.1 along the ^/-direction, (b) The contour plot 
of the three-level s' — marginals. It is very similar to (a) 
indicating the three- level system's physical dynamics may be 
considered as being primarily composed of a two-level system 
plus an inactive two-level system. However, the shifts intro- 
duced in (a) and slight differences between (a) and (b) are an 
indication of some interference between these two pictures. 



function of s f . Surprisingly, we find the same correlation 
at s = 0. This result is surprising because, although no 
inactive phase exists for y s , it is still a valid dynamical or- 
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FIG. 13. Plots of the quantum-jump activity Av, f° r different 
quadrature biases s in the three-level system. The jump ac- 
tivity grows as we bias the X-quadrature regardless of sign of 
s. However the photon count exhibits a dynamical crossover 
in the rare y-trajectories about s — 0, and shows that the 
sign of s determines the behaviour of k s /. Furthermore the 
photon variance in these quadratures also indicates the inac- 
tive/active phase crossover, with a peak at s = correlating 
with the dynamical phase transition at this point 




FIG. 14. Plots of the activity y s =o as a function of photon 
bias s' in the three-level system. y s =o exhibits a crossover at 
s f — 0, confirming that while k s i varies with s in a highly 
correlated way, y s is similarly correlated with s f . 



der parameter due to its connection to the quantum-jump 
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activity: for y s =o^ the region where s' > corresponds to 
the active quantum-jump phase and s' < corresponds 
to the inactive phase. This connection between y s= o and 
k S f also exists in the two-level system studied in Sec. 
and may be a trait present in many systems. 

To conclude, the rinding that quadrature activities may 
encode the same dynamical statistics as the photon activ- 
ity is quite surprising, given the non-commuting nature 
of the operators. We next ask whether there are sys- 
tems where a dynamical crossover is only visible through 
the quadrature activity. We will discuss this question in 
Sees. Eland ED 



V. TWO COUPLED TWO-LEVEL SYSTEMS 

We extend our study to a pair of coupled two-level 
systems driven by lasers with different polarisations, de- 
picted in Fig. [T5j In the previous two sections, we have 
identified that crossovers in quadrature activity mirror 
crossovers in the number of emitted photons, showing an 
equivalence between quadrature activity and photon ac- 
tivity as dynamical order parameters. Here we show that 
crossovers can occur in the quadrature activity and not 
in the photon activity, emphasising the use of quadrature 
activity as a dynamical order parameter in its own right. 




FIG. 15. A schematic diagram of two weakly- coupled two- 
level systems driven by lasers of opposite circular polarisation. 
The coupling between the ground states is neglected in our 
case. 

We consider two weakly coupled two-level systems as 
shown in Fig. [15l There are now two Lindblad opera- 
tors associated with each two-level system: Li = yJUci 
and Li i = \fRcn, where these operators are the ladder 
operators of each two-level system, I and II respectively. 
These two-level systems have identical Rabi frequencies 
and decay rates. As previously, we choose n = 40, but 
distinguish the two subsystems by the polarisation of the 
driving laser. Subsystem I is driven by a a x — a y po- 
larized laser, while subsystem II is driven by a <j x + a y 
polarized laser. The Hamiltonian is 



H 



ft(ci + c 
A(|0> <2| 



- ici + 
|2> (0|) 



-II + C H + lC H 



(25) 



This model allows quantum coherence between the two- 
level systems to be preserved. We choose the coupling A 



<C O, with A = Q/10 in the presented results. The weak 
coupling along with the similar physical properties of two 
subsystems is such that dynamical crossovers in quantum 
jump trajectories of this system can occur. However the 
different polarizations of the lasers driving each subsys- 
tem affords the possibility of a transition in quadrature 
trajectories. 




FIG. 16. Plots of the X-quadrature LD function as well as 
the activity and dynamical variance, for the four-level system 
depicted in Fig. [15] Note that the dynamical variance has 
been rescaled to be a factor of 5 smaller. 

While the statistics of the F-quadrature is featureless, 
Fig. [T6l shows that the X-quadrature displays a crossover 
around s = 0. This crossover is smeared out if we increase 
A and becomes sharper if we decrease A. The s < re- 
gion has a positive X-activity; this is due to light emitted 
from subsystem I. When s > 0, the light is emitted from 
subsystem II and the X-activity is negative. We expect 
that the s = behaviour is a combination of these two 
distinct dynamical phases and, to demonstrate this is the 
case, we construct time- independent probability distribu- 
tions, or phase portraits, for each dynamical regime. 

We wish to construct the marginal distributions 
e -</>O a ) for rare X-quadrature trajectories. We will bias 
the X-quadrature statistics using a field denoted s" and 
measure the typical X a for all a using the conjugate field 
s. To this end we apply the double biasing scheme as in 
Eq. ([T5]) where instead of biasing with the jump oper- 
ator and measuring X a , or vice- versa, we now bias the 
X-quadrature. This leads to an equivalent doubly-biased 
master equation 

r" (t) = C +p''"^ cos (a) 

^(e-V^eV'cl) 

-T,^p ss "+p ss "4) 



SWK ■ 



S" 2 ss>> S 2 s 

P H P 

8 F 8 F 



(26) 



where the sum is over the two subsystems. Using s" we 
bias the system to rare X-quadrature trajectories and 
we look at derivatives of the LD function with respect 
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(a) s"= -0.1 



(b) s"= 




-2-1.5-1-0.5 0.5 1 1.5 2 
X-Activity 



-2-1.5-1 -0.5 0.5 1 1.5 2 
X-Activity 



(c) s"= +0.1 



1.5 - 



f 0.5 
t5 



-1.5 




-0.5 - 1 



-2-1.5-1-0.5 0.5 1 1.5 2 
X-Activity 



FIG. 17. Phase-space portraits for the two coupled two- level systems. The X- quadrature bias in (a), (b) and (c) is s" = —0.1, 
and +0.1 respectively. In (a) and (c) the probability distributions are more concentrated about x > and x < respectively, 
whereas at s" — they are even functions of x s . These plots are indicative of the crossover in the X-activity at s — and 
highlight the less dramatic phase changes we may see using this order parameter as opposed to the jump activity. 



to s in the limit s -+ to measure the quadratures of 
these biased trajectories. Solving for the probability dis- 
tributions, for s" = —0.1, and +0.1, we find that the 
behaviour of typical (s = 0) is made of two distinct dy- 
namical phases, as is clear in Fig. [T71 

This crossover of emission from one subsystem to an- 
other is not observable in the jump trajectories and high- 
lights that quadrature activities are not just equivalent to 
the jump activity as a dynamical order parameter. These 
quadrature activities may reveal extra dynamical phases 
which are not visible simply by counting photons. 



VI. MICROMASER 

The final example we present is a many-body problem 
consisting of a set of two-level atoms interactin g w ith a 
single cavity mode, called a micromaser fl9M2ll l39j . 
We begin with a brief description of this problem be- 
fore providing a mean-field treatment to determine the 
LD function. We then proceed to determine the full LD 
function with exact numerical diagonalisation and pro- 
duce quadrature activity phase diagrams for the system. 
Finally, we examine the doubly-biased trajectory spec- 
trum of the system. We then summarise our findings, 
before moving on to present our conclusions in the final 
section of this work, Sec. IVIII 

The micromaser is a single-mode resonant cavity cou- 
pled to a thermal bath and pumped by excited two-level 
atoms which pass through the cavity. We denote the total 
number of atoms which pass through the cavity divided 
by the cavity lifetime by N ex . The steady-state cavity 
occupation distribution can change from uni-modal to 
bi-modal, depending on N ex and the atom-cavity cou- 
pling. Here, we fix N ex = 100. As we increase the cou- 
pling between the cavity mode and the atoms, we reach 
many points where the steady state cavity occupation 
number exhibits a bistability: at these points, a small in- 
crease in the coupling leads to a crossover in occupation 
number where the occupation number changes dramati- 



cally. Previously this bistability in the state of the cavity 
was shown to have an equivalent dynamical bistability in 
quantum-jump trajectories in Ref. [5]. A rich dynamical 
phase structure was found, where the number of atoms 
which have changed state while traversing the cavity was 
used to quantify the dynamical activity. 

We now turn to the model for the micromaser. Tracing 
out the atom and bath degrees of freedom, one obtains 
the superoperator, W, which contains no coherent evo- 
lution term and is of Lindblad form that is of the form 
in Eq. (fTTj) with H = and four pairs of Lindblad op- 
erators. There are two associated with the atom-cavity 
interaction 




and two result from the cavity-bath interaction, 

L 3 y/v + la (29) 
L A = yW • (30) 

Here <p is the accumulated Rabi frequency which encodes 
the atom-cavity interaction, a) and a are the cavity rais- 
ing and lowering operators respectively. Choosing a zero- 
temperature bath, we will set the thermal occupation 
number of the bath v = 0. 

We study the quadratures of the light emitted into to 
the bath and therefore deform the superoperator by the 
non-equilibrium quadrature field to 

2 

W s (p s ) = W (p s ) - | (e- ia L 3 p s + e ia L\p s ) + 8 -pP. 

(31) 

For this system, as the normal dynamics of the system 
are purely Lindbladian in nature, the LD function for the 
quadratures is identical for all a. Therefore in the follow- 
ing discussion we will focus solely on the X-quadrature 
(where a = 0) from here on. 
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A. Mean-Field Approximation 
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yja^a + 1 
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+ e |s| Vata + 1 a . 



(33) 



A variational approach to calculating the largest eigen- 
value of W s in Eq. ([33]) can be constructed using a co- 
herent state ansatz. This amounts to setting a = e %1 \fn 
and at = e~ %1 \fn and solving the corresponding Euler- 
Lagrange equations, dW s /dj = and dW s /dn = 0. 
Solving the first equation one obtains, 



FIG. 18. Quadrature activity in the micromaser from mean- 
field theory. We find multiple first order transition lines in the 
activity either side of s = 0. Also we note some bending of 
these transition lines occurs as they approach s — 0. The form 
of the mean-field theory is very similar to one obtained for the 
"atom" counting case of Ref. [5] , highlighting the connection 
between transitions in quadrature and jump-activities. 

Although a full analytic solution for the 8 = steady 
state density operator is well known 0, [39| , away from 
s = it is not easily generalised. In a study of quantum- 
jump trajectories of this system [5] progress was made by 
assuming the eigenmatrix of an ^-dependent superoper- 
ator was diagonal in at a. This approach demonstrated 
that the cavity pump rate controlled the properties of 
the coexistence line at s' = and the findings agreed 
well with exact diagonalisation results. Mulitple first or- 
der transition lines were observed in the photon active 
region where s f < 0, and a single transition was observed 
in the photon inactive regime where s' > 0, along with 
a critical point at s' ~ and </> ~ 0.1. In this work, the 
effects of quadrature biasing with s using W s in Eq. [31] 
leads to eigenmatrices which are not diagonal in at a for 
general s. However close to s = we approximate the 

non-diagonal term, — | ^Lzp s + L\p s ^, with a diagonal 

one given by 



(ap s 



P s a r ) 



, + («-I)"*(« , -|)-T' 



-ap a 



-ap s a ] + e |s| ap s a f . 



(32) 



Here we have assumed that a and at are approximately 
y/n, where n is the cavity occupation number. Thus for 
s <C 1 we may approximate the non-diagonal piece by 
an exponential and the appearance of \s\ reflects the fact 
that 9x(s) must be an even function of s. 

Through this crude approximation we have restricted 
our analysis to density operators diagonal in the number 
basis, and so the generalized quantum master equation 
reduces to an operator 



fa \e-\ s \N ex 



sin 2 (<j)y/n + l) > 



1/2 



71+ 1 



V n+ 1 



-1/2 



(34) 



Substituting these into W s ([33]) one obtains a variational 
"free energy", F s (n), whose minimum with respect to n 
yields an estimate for the LD function, 



0x0 s ) ~ -min n F s (n) . 



(35) 



This minimization process was performed numerically 
and the result is displayed in Fig. [18] Minimization re- 
veals that multiple first-order transitions occur in both 
the cavity occupation number and quadrature activity. 
These transition lines begin to bend as they approach 
s = 0. The first transition line ends at s « 0, ^ « 0.1. 
This point is the critical point discovered in Ref. 
which controls the photon number dynamics. Here it 
is much more masked, even within this crude diagonal 
approximation, than in the study of jump-trajectories of 
the atoms presented in Ref. [Ej]. Within this diagonal 
approximation, Eq. ([33]) . there are limitations and fur- 
ther implicit assumptions: the approximation becomes 
less accurate at larger (j) where non-linearities in the op- 
erator Eq. ([33]) are more prominent. There is also an 
implicit assumption of normal ordering and that the av- 
erages of products of the raising and lowering operators 
may be replaced by products of their averages. Despite 
ignoring the off-diagonal behaviour W s , this free energy 
predicts multiple first order transitions in the quadrature 
activity and, moreover, that these transitions are due to 
transitions in the occupation number, n, of the cavity. 



B. Full numerical diagonalization 

We now examine the generalized master equation in 
order to determine the exact form of the LD function. 
The generalized master equation may be expressed in 
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FIG. 19. Quadrature activity phase diagrams in quantum-jump biased systems with s' — —0.005, and +0.005 in plots (a), 
(b) and (c) respectively. In all cases the activity shows first order transition lines as we change 0. Comparing (b) with the 
mean field theory, agreement exists up to ~ 0.7, however there is significant bending of these lines as we approach the normal 
system dynamics. The degree of this bending is linked with the jump activity: in the active phase (a), the transition lines bend 
less approaching s = 0, while in the inactive jump phase, (c), the bending is more prominent. 



(a) (b) (c) 




s s s 



FIG. 20. Cavity occupation number for double counting cases with jump bias s' = —0.005, and +0.005 in plots (a), (b) and 
(c) respectively. In the jump-active case, (a), the reduced bending of the transition lines compared with (b) correlates with an 
increased cavity occupation which remains larger than the typical value away from s = 0. Conversely the increased bending in 
(c) compared with in (b) is associated with larger regions away from s — attaining typical (s = 0) cavity occupation numbers 
rapidly as we increase <j). 



terms of a matrix and diagonalized numerically. The 
method we employ is similar to that in Ref. @ , and out- 
lined in Appendix O The full phase diagram is shown in 
Figure fT9T b), where we see that the critical point previ- 
ously discussed manifests itself as the only point where 
the first order transition lines accumulate at s = 0. Be- 
yond (j) ~ 0.1 the transition lines begin to bend in such 
a manner as not to accumulate at the s = line. It 
is clear however that the mean-field theory in Fig. [18] 
predicts values of <p where the transitions occur up to 
4> ~ 0.7. These transitions correlate with the cavity- 
occupation plot shown in Fig. l20T b), demonstrating that 
these first order transitions in quadrature activity are 
also associated with the bistability of the cavity photon 
number [5|. 

We now examine the statistics of double biasing, using 



the field s' to place the system in higher or lower photon 
number state and then examining the typical quadra- 
ture realisations. Notice that the transition lines shift 
such that they either accumulate more at s = or bend 
away further before reaching s = when placing the sys- 
tem in a high (V < 0) or low (s' > 0) cavity number 
respectively. Placing the system in a fc-active phase re- 
moves lower-cavity- number regions near s = 0, visible in 
Fig. l20f b), causing the quadrature activity lines to meet 
at s = 0. Conversely, biasing the dynamics towards k- 
inactive trajectories ensures that the cavity occupation 
is small unless a large s-bias is applied. This can be seen 
in the formation of a band around s = where both 
the quadrature activity (Fig. [19]) and cavity occupation 
(Fig. [20]) are small. 

Although sharp crossovers occur in the quadrature ac- 
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tivity, corresponding to crossovers in the cavity occupa- 
tion, it is not clear that the s = dynamics are made up 
of two distinct phases. These phases are only apparent by 
examining k. Furthermore, even though the critical point 
exists in both the k and x phase diagrams, it is much eas- 
ier to discern in the former. To conclude, we find that 
quadrature measurements in micromaser can be used to 
identify dynamical phases, but in sharp contrast to the 
four-level system in [VJ the study of quadrature trajec- 
tories does not give significant insights into the system 
dynamics beyond those available by counting photons. 



VII. CONCLUSIONS 

In this paper we developed an s-ensemble approach 
to gain insight into the behaviour of light emitted from 
several quantum-optical systems. Previous studies used 
the 8-ensemble to describe the dynamics of quantum- 
jump trajectories and understand the statistics of pho- 
ton emissions in open quantum systems. We reformu- 
lated the approach in terms of generating functions and 
characteristic operators, treating the reservoir stochas- 
tically. Using this approach we examined the statistics 
of the quadratures of light emitted from various systems 
and re-constructed appropriate phase space portraits and 
Wigner functions which describe the state of the emit- 
ted light. We extended this scheme to examine the light 
not only emitted during rare quantum-jump trajectories 
where, for example, the number of emitted photons is 
larger or smaller than for average dynamics. We also ex- 
amined the quantum-jump statistics of the quadrature- 
biased trajectories, ascertaining the observable photon 
count when the quadratures of emitted light depart sig- 
nificantly from their dynamical averages. We used these 
results to understand whether different dynamical order 
parameters are correlated, and to select the most appro- 
priate s-field and order parameter with which to charac- 
terise dynamical phases. 

Our results show that even simple quantum-optical 
models such as two- and three-level systems show in- 
teresting features in the space of quadrature trajecto- 
ries. For example, for the case of the blinking three-level 
system, quadrature phase-space portraits indicate that 
emitted light may be identified with a mixture of the light 
from two simple active and inactive subsystems, which is 
consistent with the dynamical coexistence description in 
terms of number of emitted photons described in Ref. Q . 
For the case of the two-level system, our analysis here 
shows that the 'special point' in parameter space found 
in Ref. [4[ is also special for another reason: it is the 
point where one may measure the F-quadrature activ- 
ity y s in place of the quantum-jump activity k s > of the 
physical system, as they are identical in magnitude: for 
this system y s represents an alternative dynamical order 
parameter, as biasing trajectories towards larger y s leads 
to smaller £v=o, an d vice versa. 

Examination of a pair of coupled two-level systems 



showed that, in some circumstances, the quadrature ac- 
tivity reveals a phase structure which is completely unob- 
servable through photon count statistics. This highlights 
the usefulness of studying quadrature trajectories, and 
illustrates that greater inference about an open system 
can be obtained by such measurements. However, it also 
highlights that the choice of dynamical order parame- 
ter is highly system dependent. We further illustrated 
the difficulty in the choice of order parameter by study- 
ing quadratures in the micromaser. The phase diagram 
for the quadrature dynamics is rich and displays many 
features of the quantum-jump activity phase diagram 
studied in Ref. [5j, such as numerous first order tran- 
sition lines. However, the lack of direct control over the 
cavity-state bistability means one cannot simply describe 
the s = dynamics in terms of two distinct quadrature 
phases. 

Finally, we emphasize the significance of our results 
in light of homodyne detection schemes, which allow X- 
quadrature trajectories to be measured. Furthermore, we 
note that the F-quadrature statistics discussed are acces- 
sible to experimental measurement with a suitable mod- 
ification to the driving laser polarisation. We also add 
that quantum tomographic techniques allow reconstruc- 
tion of the Wigner functions of such systems. We further 
note that the methodology presented in this work can be 
extended to other more complex measurement schemes 
with other unravellings of the quantum dynamics for dif- 
ferent choices of bath measurement operators with both 
discrete and continuous spectra. 
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Appendix A: Derivation of Quadrature s-Ensemble 
Master Equation 

Using a weakly coupled Markovian reservoir one may 
consider that the ladder lowering operator of the reservoir 
acts as a driving field for the time-evolution operator of 
the total system, and so the effects of the reservoir ladder 
operator fields permit a description in terms of an input- 
output formalism In general, such a formalism works 
on the foundation that outside the range of interaction, 
the field which interacts with the system is the sum of 
the input field, that is the field just prior to interacting 
with the system, and the output field, the state of the 
field just after interaction with the system. The reservoir 
ladder operators, b(t) and b^(t f ), are the fields interacting 
with the system, the time t is the initial time when the 
fields interact with the system and, due to our choice 
of a reservoir which admits a Fock space representation, 
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one may use a stochastic description of the system-plus- 
bath dynamics. Although the input-output formalism 
is intuitive, especially if one is working with quantum 
optics where the input/output fields are manifested as 
light beams, the stochastic description [3| of the system 
allows for a concise mathematical representation where 
the effects of the reservoir on the system are reduced to 
an Ito increments which obey a quantum Ito calculus. 

Before we proceed with the derivations of the general- 
ized master equation for biased quadrature ensembles, 
it is necessary to outline the calculus of the reservoir 
increments, as well as the quantum-jump dK process. 
Working with an unsqueezed vacuum, the only non-zero 
combination of dB(t) and dB^(t) is: 

dB(t)dB ] (t) = dt. (Al) 

With this Ito table we can deduce the appropriate table 
for the jump trajectory process dK, 

dK{t)dK(t) = dK(t) 
dB(t)dK(t) = dB(t) 
dK(t)dB\t) =dB\t). (A2) 

Now, using the definition of p s ([T3j), we examine the in- 
crement of the s-biased density matrix 

d[p s ] = Tr Res (d[V£ (t)]p + V$ (t) d[p] + d[V£ (t)]d\p]) . 

(A3) 

The first two terms appear in standard calculus while the 
final term appears due to the stochastic nature of these 
increments. We have already defined dp in Eq. (fT0]h and 
so we need to examine the increment of the characteristic 
operator. Expanding (£), defined in ([T5]h and using 
the calculus set out in Eq. (jAl|) one finds 

d[VZ„ (t)} = (t) (jdt - S -dX<^j . (A4) 

We may now evaluate the expansion on the right-hand 
side of (|A3|h term by term, to obtain 

Tr Res (d\V$ a (t)]p) = jp s dt, 

(t)d[p]) = C(p s )dt, 
TrRes (d[V£ a (t)]d[p}) = - J2 \ {e~ la L iP s + e ia p s L\) dt, 

(A5) 

where the Li are the Lindblad operators of the system 
coupled to the environment. We can now see the incre- 
ment equation (jA3|) leads to Eq. (fT7|) . 

This procedure may be extended straightforwardly to 
the doubly-biased cases. In the case of the biasing pho- 
ton trajectories and examining the typical quadrature 
statistics (Eq. ([18]) ) one introduces the jump-trajectory 
characteristic operator V£ (t) as in Eq. (|A1|) and defines 
its increment 



«' (*)] = V S K (t) (e" 5 ' - l) dK. (A6) 

Now we proceed as before using the Ito calculus set out 
in Eqns. ([AT]) and (lA2l) to arrive at the result Eq. (fl8j) . 
Finally we note that, although introduced for a specific 
set-up, this formalism is not restricted to pure states and 
allows a freedom of choice for the reservoir provided we 
stay within the Markovian approximation. For exam- 
ple, it may be readily generalized to a squeezed thermal 
reservoir. 



Appendix B: Marginal Distributions, the Wigner 
Distribution and the Inverse Radon Transform 

In the analysis of quadrature statistics, we consider 
the set of marginal distributions e~^ x ^\ associated with 
the general quadratures defined in Eq. (|7)), where a 
varies over the range [0,7r]. We also employ the quasi- 
probability distribution of the light, known as the Wigner 
distribution [J-Q- W e use these distributions to study 
the light in the reservoir coupled to the system. 

Taking the density matrix of the reservoir, p res , the 
Wigner function in the Schrodinger representation is 

W„{W) = ±J d 2 6e- M ' +s ^Tr ies (p re J s » 1 ~ s ' b )) . 

(Bl) 

The variable d is identified with the quadratures through 
the relations 

X = Re (#) (B2) 
Y = Im (0) . (B3) 

As direct evaluation of this distribution is difficult we will 
use the fact that W p is related to the marginals e"^^ 
by the Radon transform [io|, via 

e -^ a ) = n [ Wp ] (B4) 

/+oo 
Wp (X cos a — Fsina, Fcosa+Xsina) dY 
-oo 

and we employ its numerical inverse to the marginals 

to obtain W p . This technique is frequently used in quan- 
tum tomography Q. The Wigner distribution allows 
classification of the light into quantum and classical com- 
ponents based on the sign of the distribution: negative 
regions are indicative of quantum effects such as entan- 
glement. 
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Appendix C: Full diagonalization of the s-biased 
generalised master equation for the micromaser 

To determine the LD function of the micromaser nume- 
rially, we numerically diagonalized the generalized master 
operator W s in matrix form, as in Ref. [7] . It is necessary 
to truncate the basis of the system. The basis of num- 
ber states \n) is suitable and, for the N ex value of 100 
studied in this work, we restrict the maximum photon 
number to n = 150. The form of W s introduces co- 
herences between different eigenstates of the cavity and 
so we must ensure our basis allows these coherences to 
be preserved. However, since in matrix form the rep- 



resentation of the superoperator is an n 2 x n 2 matrix, 
we investigated whether further truncation was possible. 
In practice it is possible to truncate the basis in such a 
manner that only coherences between number states with 
occupation number differing by m (with m < n) are pre- 
served. The value of m may be tested numerically to 
ensure the results are not sensitive to this truncation. In 
the results presented in this work, the values of n = 150 
and m = 15 were found to be sufficient for the diago- 
nalization process. The largest real eigenvalue extracted 
from this matrix is identified as (s). The matrices were 
diagonalized using an Arnoldi iteration scheme [41]. The 
activity was calculated by taking numerical derivatives 
of the resulting LD function. 
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